In [1]:
pwd
Out[1]:
'C:\\Users\\alpha\\Documents\\TestingGeopandas'
In [2]:
cd C:\Users\alpha\Documents\work\burned pixels\ecoregions_countries_climates_continents
C:\Users\alpha\Documents\work\burned pixels\ecoregions_countries_climates_continents
In [3]:
import geopandas as gpd
import matplotlib.pyplot as plt
import pickle
import pandas as pd
import numpy as np
import calendar

import matplotlib.cm as cm
import matplotlib.colors as cls
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
In [4]:
#latex
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'custom'
matplotlib.rcParams['mathtext.rm'] = 'Times New Roman'
matplotlib.rcParams['mathtext.it'] = 'Times New Roman:italic'
matplotlib.rcParams['mathtext.bf'] = 'Times New Roman:bold'
In [5]:
from scipy.stats import boxcox, probplot
In [6]:
bios = gpd.read_file('biomes_continents')
In [7]:
bios.columns
Out[7]:
Index(['OBJECTID', 'BIOME_NUM', 'CONTINENT', 'is_north', 'SUM_2001_J',
       'SUM_2001_F', 'SUM_2001_M', 'SUM_2001_A', 'SUM_2001_1', 'SUM_2001_2',
       ...
       'SUM_2017_3', 'SUM_2017_4', 'SUM_2017_S', 'SUM_2017_O', 'SUM_2017_N',
       'SUM_2017_D', 'SUM_Area', 'Shape_Leng', 'Shape_Area', 'geometry'],
      dtype='object', length=212)
In [9]:
area = bios['SUM_Area']*1e-6
In [12]:
data_columns =bios.columns[4:4+17*12]
In [13]:
inter_cols = {}
for yr in range(17):
    inter_cols[yr] = data_columns[12*yr:12*yr +12]
In [14]:
intra_cols = {}
for mo in range(12):
    intra_cols[mo] = [data_columns[mo + yr *12  ] for yr in range(17)]
In [15]:
for yr in range(17):
    inter = bios[  inter_cols[yr]].sum(axis=1)/area
    label = 'inter_' + str(yr +2001)
    bios= bios.assign(label=inter)
    bios = bios.rename(columns={'label': label})
In [16]:
inter_columns = bios.columns[-17:]
In [17]:
inter_max = bios[inter_columns].max().max()
In [18]:
for mo in range(12):
    intra = bios[ intra_cols[mo]].sum(axis=1)/(17*area)
    label = 'intra_{:02}'.format(mo)
    bios = bios.assign(label=intra)
    bios = bios.rename(columns={'label': label})
In [19]:
intra_columns = bios.columns[-12:]
In [20]:
intra_max = bios[intra_columns].max().max()
In [21]:
def draw_map(df,column,title=None,vmin=0,vmax=1,name='name',dpi=100):
    fig, ax = plt.subplots(figsize=(12, 8))

    ax.set_aspect('equal')
    ax = df.plot(ax =ax, column=column, cmap='OrRd',vmin=vmin, vmax=vmax)
    ax.set_axis_off()
    
    norm = cls.Normalize(vmin, vmax)
    cmmapable = cm.ScalarMappable(norm, 'OrRd')
    
    cmmapable._A = []
    cbaxes = inset_axes(ax, width="80%", height="1%", loc=3)
    cb = plt.colorbar(cmmapable, cax=cbaxes,orientation="horizontal") 

    cb.set_label(title, fontsize=20, family='Times New Roman')
    cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=15, family='Times New Roman')

    plt.show()
    fig.savefig(name, dpi=dpi, bbox_inches='tight')

Inter-anual

In [22]:
lambdas = []
for yr in range(17):
    column= 'inter_{}'.format(2001 + yr)
    test_column, lbd = boxcox(bios[column] + 1e-25)
    lambdas.append(lbd)
lmbda=  sum(lambdas)/len(lambdas)
In [23]:
lambdas
Out[23]:
[0.09794875021807274,
 0.10709764548918331,
 0.11000933794034778,
 0.11648864557704566,
 0.11049426552203133,
 0.1136208500903649,
 0.11947258901139729,
 0.10740846861730337,
 0.10432852879128167,
 0.11260004417237568,
 0.08974439411896297,
 0.0999864366237989,
 0.10997392807044461,
 0.10670316644737485,
 0.10572148632865014,
 0.11486313745614565,
 0.11754796825230394]
In [24]:
probplot(bios.inter_2017, plot=plt.gca());
In [25]:
probplot(test_column, plot=plt.gca());
In [26]:
maxs = []
mins = []
for yr in range(17):
    column= 'inter_{}'.format(2001 + yr)
    test_column = boxcox(bios[column] + 1e-25,lmbda=lmbda)
    maxs.append(test_column.max())
    mins.append(test_column.min())
In [27]:
box_min = min(mins)
box_delta =max(maxs) - box_min
In [28]:
for yr in range(17):
    column = 'inter_{}'.format(yr +2001)
    test_column = (boxcox(bios[column] + 1e-25,lmbda=lmbda) - box_min)/box_delta
    bios = bios.assign(test_column =test_column)
    yr_name = str(yr +2001)
    #"Burnt area pixels (#/km2) by biome in 2001
    title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome parts in ' +yr_name
    #title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
    #title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
    name = 'inter_biomes_parts_index_{}'.format(yr +2001)
    draw_map(bios,'test_column',title=title,name=name)
C:\ProgramData\Anaconda3\envs\geopandas\lib\site-packages\matplotlib\font_manager.py:1339: UserWarning: findfont: Font family ['cursive'] not found. Falling back to DejaVu Sans
  (prop.get_family(), self.defaultFamily[fontext]))

Intra-anual

In [29]:
lambdas = []
for mo in range(12):
    column= 'intra_{:02}'.format(mo)
    test_column, lbd = boxcox(bios[column] + 1e-25)
    lambdas.append(lbd)
lmbda=  sum(lambdas)/len(lambdas)
In [30]:
lambdas
Out[30]:
[0.08468620438387275,
 0.08432306246379806,
 0.0982259003412448,
 0.08332125392466332,
 0.0892014445845387,
 0.09068062588895898,
 0.10526759762768141,
 0.10205923524100768,
 0.10482478599652979,
 0.099680805112141,
 0.08900407966603739,
 0.07327253953187773]
In [31]:
probplot(bios.intra_11, plot=plt.gca());
In [32]:
probplot(test_column, plot=plt.gca());
In [33]:
maxs = []
mins = []
for mo in range(12):
    column= column= 'intra_{:02}'.format(mo)
    test_column = boxcox(bios[column] + 1e-25,lmbda=lmbda)
    maxs.append(test_column.max())
    mins.append(test_column.min())

box_min = min(mins)
box_delta =max(maxs) - box_min
In [34]:
for mo in range(12):
    column = 'intra_{:02}'.format(mo)
    test_column = (boxcox(bios[column] + 1e-25,lmbda=lmbda) - box_min)/box_delta
    bios= bios.assign(test_column =test_column)
    mo_name = calendar.month_name[mo + 1]
    #"Burnt area pixels (#/km2) by biome in Janury (2001-2017)"
    title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome parts in ' +mo_name + ' (2001-2017)'
    #title = r' Burned Area Pixels/$\mathrm{Km}^{2}$ by climate 2001-2017 ' + mo_name
    #title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
    #title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
    name = 'intra_biomes_parts_index_{:02}'.format(mo)
    draw_map(bios,'test_column',title=title,name=name)

Total

In [35]:
total = bios[data_columns].sum(axis=1)/17

total = total/area

bios= bios.assign(total=total)
In [38]:
column = 'total'
test_column, lbd = boxcox(bios[column] + 1e-25)
test_column = (test_column - test_column.min())/(test_column.max() - test_column.min())
bios = bios.assign(test_column =test_column)
title = r'Burned area index (#/$\mathrm{km}^{2}$) by biome parts in (2001-2017)' 
#title = r'Burned Area Pixels/$\mathrm{Km}^{2}$ by Climate 2001-2017 {} '.format(calendar.month_name[mo + 1])
#title =(r'ABC123 vs $\mathrm{ABC123}^{123}$')
name = 'total_biomes_parts_index'
draw_map(bios,'test_column',title=title,name=name)